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ABSTRACT 

We present starburst models for far- infrared/sub-millimeter/millimeter 
(FIR/sub-mm/mm) line emission of molecular and atomic gas in an evolving 
starburst region, which is treated as an ensemble of non-interacting hot bub- 
bles which drive spherical shells of swept-up gas into a surrounding uniform gas 
medium. These bubbles and shells are driven by stellar winds and supernovae 
within massive star clusters formed during an instantaneous starburst. The un- 
derlying stellar radiation from the evolving clusters affects the properties and 
structure of photodissociation regions (PDRs) in the shells, and hence the spec- 
tral energy distributions (SEDs) of the molecular and atomic line emission from 
these swept-up shells and the associated parent giant molecular clouds (GMCs) 
contains a signature of the stage of evolution of the starburst. 

The physical and chemical properties of the shells and their structure are 
computed using a a simple well known similarity solution for the shell expansion, 
a stellar population synthesis code, and a time-dependent PDR chemistry model. 
The SEDs for several molecular and atomic lines (^^CO and its isotope ^^CO, 
HON, HCO+, C, O, and C^) are computed using a non-local thermodynamic 
equilibrium (non-LTE) line radiative transfer model. 

By comparing our models with the available observed data of nearby infrared 
bright galaxies, especially M 82, we constrain the models and in the case of M 82, 
we provide estimates for the ages (5-6 Myr, 10 Myr) of recent starburst activity. 
We also derive a total H2 gas mass of ~ 2 - 3.4 x 10^ Mq for the observed regions 
of the central 1 kpc starburst disk of M 82. 



Subject headings: ISM: clouds - ISM: evolution - galaxies: ISM - galaxies: star- 
burst - galaxies: star clusters - radio lines: ISM 
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Introduction 



Knowledge of the physical properties and evolution of the gas and dust content in 
the interstellar medium (ISM) of starburst galaxies is essential for understanding the cause 
and temporal evolution of star-forming activity. In particular, studies of such galaxies in the 
nearby universe are essential as a step in understanding the role of the starburst phenomenon 
in the cosmic evolution of galaxies. To constrain theories of how the ISM evolves, one 
needs to investigate both individual galaxies and large statistical samples of data at multiple 
wavelengths. Especially, with the available data for the dust component, studying the gas 
in the co-space ISM becomes more interesting and important. 

The number of multiwavelength observations of starburst galaxies throughout the cosmic- 
scale has increased dramatically due to the significant improvement in the sensitivity and 
resolution of telescopes. These observations provide an essential basis for starburst modeling, 
and such modeling provides systematic predictions of the properties of the ISM in idealized 
models of starburst galaxies for comparison with these observations. For example, M 82, 
an irregular starburst galaxy (10) in the nearby universe has an inclination of 81°. The in- 
frared luminosity of M 82 arises mostly from the central 1 kpc region, which has a stellar bar 



struc t ure and a comp lex system of clumps and filaments (ILynds fc Sandagelll963l : iBeck et al. 



19781 : ILo et al.l 119871 ) . The formation mechanism of this complex system and the evolution- 
ary scheme in M 82 remain under deba te (e.g. IVisvanatharu 119741 ICarlstrom fc Kronberg 
199ll : IShen &: Ldll996l : IWills et al.ll2000l ). Most recently, we conducted an ideal case study 
of an expanding shell model in M 82, and suggested that the circumnuclear rings seen in this 
galaxy may possibly be a consequence of the evolution of swept-up gas caused by starbursts 
that occurred in the center ~ 100 Myr ago (Yao et al 2006, hereafter Paper I). 

Previously, the age estimates of the starburst in M 82 have been presented by many 
authors. Yun et al. (1993) compared the disk HI with optical maps, and found a large 
amount of gas being channeled into the core of the galaxy over the last 200 Myr due to the 
tidal encounter with its large spiral neighbor galaxy M 81 . Ages derived from super star 
clust ers (SSCs) in optical and near-IR ima ges are ~ 50 Myr (Ide Griis et al.ll200ll), ~ 30 - 100 



Myr (IRieke et al 



1993 



2003: iSmith et al. 



12006 



Barker et al 



Melo et al 



20081) . ~ 4 - 6, and 10 - 30 Myr fiForster-Schreiber et al. 
2OO5I ) . The corresponding estimates of the average star 



formation (SF) rate over the 200 Myr period is roughly about 10 Mq yr~^. The importance 
of optical and near-IR spectroscopy in studies of dusty star-forming galaxies has long been 
recognized. But studies of young stellar populations at these wavelengths remain difficult. 
The age determinations are affected by residual effects due to the age-metallicity degeneracy, 
and age-IMF degeneracy. In addition, the completeness of the sampled stellar population is 
affected by the unavoidable effects of extinction in the optical and near-IR. 
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The star formation history of M 82 ha s also been studied us i ng mid-infrared and 



far-infrared spectroscopy (jCqlbert et al.l Il999l : lEfstathiou et ahl I2OOOI : iGaUiano et ahl 12003 



Siebenmorgen fc Kriigell 120071 ) . At these wavelengths, the fine structure line emissions are 



relatively insensitive to extinction, and hence can provide a unique probe of age and star for- 
mation history in an infrared-bright, dust-obscured galaxy like M 82. Colbert et al. (1999) 
obtained a burst age of 3 - 5 Myr for the central 1 kpc (65 - 85 ) region using an instanta- 
neous starburst model and a steady-state PDR model. However, their single burst model is 
dominated by the brightest and most recently formed stars (the hot spots seen in mid-IR). 
Efstathiou et al. (2000) presented an evolving starburst model for dusty media using state of 
the art codes for calculating the radiative transfer in dust shells, and incorporating a model 
for the composition and size distribution of grains in the ISM. Their study concluded that 
it is possible to relate the observed infrared spectrum of dust associated with a starburst to 
its age and its star formation history by following the evolution of an ensemble of GMCs 
of identical mass induced by massive star formation in their centers. They show that the 
burst age for the central 500 pc region of M 82 is between 10 and 30 Myr using a model with 
two instantaneous bursts. Their derived a ges are supported by near-IR spectroscopy and 
high-resolution imaging of stellar clusters (ISatyapal et al.l 119971 ). Efstathiou et al. (2000) 
also suggested that far-IR surveys may preferentially detect older starbursts than mid-IR 
studies, based on an argument concerning the evolution of the luminosity of starbursts ob- 
served at different wavelengths. It is clear that the estimation of M 82 starburst age has a 
large uncertainty. The questions we ask are 



Does the FIR/sub-mm/mm molecular and atomic line data contains a signal of the 

evolutionary phase of a starburst? 

Can the molecular sub-mm and mm lines provide an alternative tool for estimating the 

starburst ages? 



This study is motivated by the abundant evidence of giant bubbles and shells found in 
Milky Way, 30 Doradus in the Large Mag ellanic Cloud (LMC), and nearby starburst galaxies 
( jPedlar et al.l l2003l : Ide Grijs et al.l l200ll ). the success of using the dusty starburst model 
to constrain the star fo rmation history of obs erved IRAS starburst galaxies by following 
an ensemble of GMCs (lEfstathiou et al.l I2OOOI ). and the available multiple transitions in 
several molecular tracers. In a previous paper (Paper I), as a preliminary approach we 
presented a simple model of an expanding supershell surrounding a massive star cluster and 
its expansion at different evolutionary stages that predict sub-mm/mm and far-IR emission 
line intensities from several key molecular and atomic constituents in the shell (^^CO and 
its isotope ^^CO, HCN, IICO+, C, O, and C+). In this paper, we present a series of models. 
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called Evolving Starhurst Model (ESbM) for the molecular and atomic line emission from an 
evolving starburst region, which is treated as an ensemble of evolving GMCs each of which 
is centrally illuminated by a compact star cluster (SC). The GMCs in our model ensemble 
follow a power-law mass spectrum. By comparing with the available observed data on nearby 
starburst galaxies, we can constrain the models and provide better interpretations for the 
observations. The main goals of this paper are (1) to show that it is possible to model the 
FIR/sub-mm/mm line emission of molecular and atomic gas by following the evolution of a 
starburst region, as in certain infrared (IR) models; (2) to relate the observed molecular line 
properties of a starburst galaxy to its age, and hence to constrain the global star formation 
history; (3) to understand the formation mechanism of the molecular rings in M 82; (4) 
and finally, to provide useful information for the interpretation of future high resolution 
maps of molecular gas on small and large scales in starburst galaxies, in order to provide a 
deeper understanding of the structure, dynamics, and evolution of the neutral ISM and its 
relationship with active star formation. 

The plan of this paper is as follows. In § [2] we discuss briefly our model assumptions, 
initial parameters, and model procedure. In §[3] we present examples of our modeling results 
for an ensemble of GMCs/shells. In § H] we present an application of the model to derive 
the age of the starburst and molecular gas swept up by the shells in M 82, and to provide 
new insights into the nature and physical state of the ISM in its starburst region. We also 
discuss applications to the supershell surrounding supernova remnant (SNR) 41.9 + 58. In 
the conclusions, we summarize the main results of this work. 

Through this study, we hope to provide some answers to those intriguing questions 
mentioned above, and to lay a foundation for future starburst modeling for neutral gas 
media. 



2. Evolving Starburst Model For Gas Media 



The basic assumptions for our ESbM model are (1) star fo rmation occurs within the 
dense optically thick spherical cloud (GMC) (e.g. iGao et al.ll200ll ). and star formation takes 
place instantaneously, with the star cluster treated as a point source (see Paper I); (2) 
absorption of the starlight from the central cluster is produced by dust associated with the 
gas, assuming a constant gas-to-dust ratio; and (3) the gas responding to star formation in 
a starburst galaxy is treated as an ensemble of GMCs with different initial masses, each of 
which responds to massive star formation at its center. 



Our ESbM model incorporates a standard similarity model for the bubble/shell struc- 
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trans fer model for molecular and atomic lines (iRawlings fc Yatesll200ll : Ivan Zadelhoff et al. 



ture around a young star cluster (IWeaver er al.l 119771: iMcCrav fc Kafatod 119871 ). a time- 
dependent stellar population sy nthesis model (Leitherer et al. 1999[) . a fully time-dependent 



chemistry model for the PDRs JBcU et ahl bood: iRolhg et al.ll2007h . and a non-LTE radiative 



20021 ). Theoretical background for each of these four physics elements contained in our 
model has been presented in Paper I. Simulation methodology for an expanding supershell 
model has also been described in Paper I. In this paper we present an ensemble of ideal 
three-dimensional, spherical expanding shells, in order to model the line emission of neutral 
ISM in massive star-forming regions in a starburst galaxy like M 82. Our instantaneous 
starburst model does not address issues related to more complicated geometry in order to 
understand how these shells are distributed in a galaxy, how they interact, or how the gas 
becomes available for fueling the massive star formation in our model GMCs. 

In this paper, we divide the evolutionary scheme of the expanding shell/GMC ensemble 
into two phases referred to as Winds and post-SN. The Winds phase begins with the formation 
of a star cluster and an H II region inside the GMC followed by the formation of a rapidly 
expanding hot bubble produced by strong stellar winds. This phase ends when the bubble 
breaks out of its parent GMC. In this phase, the parent GMC is assumed to be stationary, 
and acts as a dense uniform ambient ISM to the expanding shell formed by gas swept up by 
the bubble. The post-SN phase starts when the most massive star in the ensemble reaches its 
main-sequence lifetime, and explodes into a supernova. The shell expansion in this phase is 
first driven by repeated supernova explosions, then changes from pressure-driven (adiabatic) 
to zero pressure {non- adiabatic) as the hot bubble cools. In the post-SN phase, the shell 
expands into a less dense uniform ambient ISM. A steady-state mechanical power and energy 
for each phase is assumed, in order to satisfy the requirements of the similarity model. Our 
non-LTE line radiative transfer method simply sums the line emission from the model shells 
and parent GMCs in the ensemble for the Winds phase, and just the model shells for the 
post-SN phase. Hence, our model is an idealistic approximation for a starburst galaxy; it 
may be considered the first step toward simulating the response of the gas environment in 
an evolving starburst region for the purpose of examining the effects of this evolution. 

A schematic diagram of the model components (PDR or shell and its parent GMC) is 
shown in Fig. [TJ In this paper, we treat the PDR and shell as one gas component. The 
line intensity/flux for the shell and GMC components are calculated using the same method. 
For the Winds phase, the integrated line intensity/flux at each time step is the sum of 
line emission in the shell and its parent cloud. For the post-SN phase, the integrated line 
intensity /flux is the emission from the shell only. Due to the incomplete knowledge of the 
structure and physical state of the ambient ISM in a starburst galaxy, we do not include the 
molecular or atomic line emission from this component in our model. The justiflcation of the 
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ambient ISM density (30 cm ^) will be discussed in § 12.2. 2[ In § |U we discuss the possible 
impact of this exclusion on the conclusions, specifically regarding the applications to M 82. 

2.1. Model Outline 

Our basic model comprises a series of non-overlapping (i.e. non-interacting) spherical 
shells expanding into a uniform gas medium. These shells are all driven by winds from star 
clusters formed during an instantaneous starburst. The interior hot bubble is produced by 
stellar wind from an underlying super star cluster, whose properties are selected as discussed 
later. The thrust of our model simulation is to compute the molecular line emission from the 
swept-up shells and the associated parent GMCs. Since the underlying stellar radiation from 
the clusters has a pronounced effect on the properties of PDR regions of the shells, and since 
these properties are therefore also affected by the radius of the shell and evolutionary stage 
of the cluster, the SED of the molecular line emission from these shells contains a signature 
of the stage of evolution of the starburst. This variation with time, predicted by our model, 
offers a way of dating the starburst, at least in principle. 

The set of our starburst models is divided into two phases, namely the Winds and 
post-SN phases. In the Winds phase, the shells propagate into their parent clouds, which are 
substantially more dense than the surrounding ISM in which they are embedded. In the post- 
S'A'^ phase, the shell breaks out of the parent cloud and expands into the uniform lower density 
ISM which pervades the entire galaxy. The same bubble/shell dynamical theory is used to 
describe the shell behavior in both phases. Since the simple similarity relations do not apply 
to a nonuniform ambient medium, we do not follow the shell expansion across the transition 
from cloud to the surrounding ISM. Instead the two phases are treated independently, and 
can be viewed as simply alternative scenarios. Continuity at the transition is only maintained 
in the mechanical luminosity of the wind (L^) and the stellar luminosity {L^,) evolution of 
the central star clusters. There is accordingly a discontinuity in the radius (Rg), velocity 
(Vg), and consequently the temperature (T^), density (n^), thickness (dg), and the mass 
of each expanding shell across the boundary between the two phases. These quantities 
asymptotically approach those of a continuous model when the mass of the ISM swept up in 
the post-SN phase becomes greater than the mass of the parent GMC in the Winds phase. 
The Winds phase thus comprises younger and denser shells than in the post-SN phase. In 
principle, it is possible to model data with starbursts occupying a large range of potential ages 
and molecular gas excitation conditions, especially if the dense molecular core component 
were also included our model. For this study, our Winds phase model results are used only 
for a comparison with the observed expanding supershell surrounding SNR 41.9 + 58 in the 
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Fig. 1. — This diagram illustrates the structural components associated with a single star 
cluster within our model. The white region is the hot cluster wind, the blue region is the 
shell of material swept up from the giant molecular cloud, represented by the orange region. 
The region exterior to the GMC is the ambient interstellar medium (ISM) with a fixed H2 
density of 30 cm~^. 
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center of M 82, which will be discussed later in § HI 

In each phase, the shell structure is computed with time as the independent variable. 
The final output dependent variables are the line fiuxes (and profiles) for several molecules 
and atoms each at a number of observed transitions, computed by a non-LTE line radiative 
transfer code applied to each shell and its parent cloud. The integrated line fiux for each shell 
(and GMC) is the sum over the emission from the entire emitting region. The total line fiux 
for the shell ensemble is then the sum of the integrated line fiuxes of all shells (and GMCs) . 
Intermediate variables which determine these fiuxes include the radius and velocity of each 
shell, its chemical structure, shell temperature and density structure, which are computed 
using a Shell Dynamics code, a stellar population synthesis code, and a time-dependent PDR 
code. These codes are described in detail in the Ph.D. Thesis of Yao (2009, hereafter Yao 
Thesis, and references therein). 

Our two-phase starburst model described above must also be characterized by a number 
of fixed parameters with adopted plausible values. These include, for example, the initial 
giant molecular cloud parameters (mass Mgmc, initial H2 density Uq, and core H2 density 
He), the star formation efficiency [r]), the star cluster related parameters (IMF, individual 
star mass m*), the initial chemical composition of the parent clouds, and the density of 
the ambient ISM. These parameters, along with others, and their numerical values, are 
discussed in detail in the subsequent sections, in association with the discussion of the PDR 
and radiative transfer codes. A brief summary of all variables and parameters are presented 
in Table [3] at the end of this section. 

Finally, a chi-square (x^) method will be used for fitting the model line spectral energy 
distribution to a set of molecular line data in order to estimate the starburst age(s) and 
observed total H2 mass in the nuclear disk of M 82. In a subseque nt paper (Yao 2 009), we 



will present applications to luminous infrared galaxies (LIRGs) (e.g. lYao et al.ll2003l ) beyond 
M 82 using our FIR/sub- mm/mm starburst model. We discuss the relationships between 
the excitation of CO molecule and SF properties of LIRG galaxies, and derive the behavior 
of the model CO-to-H2 conversion factor X in a starburst galaxy. 



2.2. Model Input Parameters 

2.2.1. Winds Phase 

In 2005 Keto et al. observed ^^C0(2-l) emission in the center of M 82 with a linear 
resolution of 17 pc at the source. They resolved ~ 300 molecular clouds with masses ranging 
from ~ 2 X 10^ to 2 x 10^ Mq. The mass spectrum of these GMCs scales as dN /dMcuc 
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oc Mj2^^°'^, similar to the Galactic one (ISanders et al.lll985l : ISolomon et al.lll987l ). Keto et 
al. also found that the mass spectrum of star clusters in M 82 follow the same power-law 
distribution, suggesting that individual molecular clouds are transformed in the starburst 
into individual star clu sters in their dense cores. Combining this result with other studies 



Weiss et all 120051 ). we assume the GMC mass distribution responsible for the stellar 



(e.g. 

outburst in our model has a power-law index of 1.5, and the mass ranges between 3.16 x 10'^ 
and 10'' Mq. About 70% of molecular gas mass in a model starburst will then be contained 
in the clouds with masses > 10^ Mq. It is also expected that much of the FIR luminosity 
due to star formation would arise from these massive clouds. To reduce the computation 
time, a discrete and arbitrary number of giant molecular clouds distributed similar to that 
discussed in Keto et al. (2005) is assumed. The masses for these discrete GMCs are 3.16 x 
10^, 10^, 3.16 X 10^, 10^ 3.16 X 10^ 10^ 3.16 x 10^ and 10^ Mq, and are hereafter denoted 
as 3M3, M4, 3M4, M5, 3M5, M6, 3M6, and M7. The total number of clouds is about ~ 
400, and the total H2 mass contained in the GMCs and shells is ~ 1.69 x 10^ Mq with a 
total star clusters mass of 4.2 x 10^ Mq. This selection is intended to provide a template 
cloud/cluster mass for scaling the model to fit the data for M 82. The best fitting molecular 
II2 gas mass and the initial star cluster mass will be determined from a fitting method, 
as described later in § HI 

The average gas densities of GM Cs in our Ga l axy and starburst galaxies a. r e in the range 



a few 10 to a few times 10^ cm ^ ( iDame et al.l Il986l : iJog &: SolomonI Il992l : IWilson et al. 



20081 ). but their cores, where most of the stars form, have much higher densities. Higher gas 
dens ities are expected in more actively star-forming galaxies in accordance with the Schmidt 
law (IKennicuttl 119981 ) . Therefore, we adopt a value of 300 cm"^ for the uniform initial H2 
gas density (i.e. no) for the M7 cloud based on the densities for the most massive clouds 
in the study of molecular cloud properties in the active spiral M 51 by Scoville & Wilson 
(2004). Since this work is proposed to be a complementary study to the dusty starburst 
models developed by Efstathiou et al. (2000), we adopt the same core density namely Uc = 
2 X 10^ cm~^ for the calculations of the number of ionizing photons (used in computing the 
Stromgren radius Rs value). The radius of this 10^ Mq cloud is 47 pc derived from the mass 
of th e cloud and the assunied density with the assurn ed effective molecular weight ^ = 2.36 
(e.g. 



Elmegreen et al.lll979l : iMcCrav fc Kafatoslll987h . 



We know that star formation takes place primarily in the dense cores of GMCs, but 
the details of the physical processes involved are not yet well understood. The efficiency 
of star formation (or the gas consumption rate) ranges from abo ut 1% in late- type spirals 
to 60% or more in active star- forming galaxies (jKennicuttl Il998l ) . In this paper we adopt 
a moderate star formation efficiency rj = 25% for our model starburst galaxies. A Salpeter 
IMF is assumed, i.e. dN/dm^ oc m~^'^^ (IMF; ISalpeterlll955l ). and the stellar mass is in the 
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range 0.1 - 120 Mq. 

We assume that the relationship between cloud mass and radius is the same as that 
derived from a CO survey for 273 giant molecular clouds in the Galactic inner disk by 
Solomon et al. (1987). From the measured relationship between the cloud size and the 
velocity line width, and the application of the virial theorem, they derived a power-law 
cloud density and mass relation, in which the mean gas density of the cloud is inversely 
proportional to the cloud size. Hence, the cloud mass is proportional to the square of the 
cloud radius (i.e. the mass surface density is a constant). From the studies of independent 
methods of determining the H2 mass, Solomon et al. (1987) also demonstrated that these 
giant molecular clouds are bound principally by self-gravity and not by external pressure 
exerted by a hot phase of the ISM. Since we assume the mass distribution of GMCs in a 
starburst galaxy is similar to that in our Galaxy, we adopted the power-law relations of mass, 
radius, and density for our model GMCs as those defined in Solomon et al. (1987). This 
mass - radius relation has also been studied for GMCs in other galaxies, for example M 51 
by Bastian et al. (2005), and has been found to be similar in form. Hence, the density and 
the radius for a model CMC having mass less than 10'' M0 can be written in the following 
forms. 



where Rgmc is the radius of the CMC with mass Mgmc that is less than 10^ M©. Table □ 
summarizes the number distribution of the GMCs and their initial physical properties. For 
CMC mass less than 3.16 x 10^ M©, the predicted number of very massive stars (> 30 M©) 
in the star cluster is below 1.0. In addition the supernova wind will not be steady as assumed 
in our model, because of the relatively small numbers of contributing stars. Initially, the 
total CMC and cluster masses are assumed to be ~ 1.86 x 10'^ M© and ~ 4.27 x 10^ M©, 
respectively. 



All shells are propelled into a less dense ambient ISM during the post-SN phase. The 
intercloud medium of the central 1 kpc region of the Galaxy has been studied by Jog & 
Solomon (1992), who find it to be mostly molecular with density between 30 and 100 cm"^. 




(1) 



(2) 



2.2.2. Post-SN Phase 



Table 1. Initial conditions of GMCs and SCs in a modeling starburst system. 



GMC (Nnorm) 


logloMcMC 




d 

ric 


Rgmc" 


logioMsc 








1 r * h 


1 r TTi, fch,\ 
logioige 




(Mo) 


(cm^'^) 


(cm^"^) 


(pc) 


(M©) 


(O.1M0 < M, < I2OM0) 


(M* > 8Mq) 


(M* > 3OM0) 


(erg s-l) 


(erg s-l) 


M7 (1.0) 


7.0 


300 


2000 


46.8 


6.4 


7.1E6 


2.2E4 


5.0E3 


42.8 


40.1 


3M6 (1.77) 


6.5 


534 


3558 


26.3 


5.9 


2.2E6 


7.0E3 


1.6E3 


42.3 


39.6 


M6 (3.1) 


6.0 


949 


6325 


14.8 


5.3 


7.1E5 


2.2E3 


5.0E2 


41.8 


39.1 


3M5 (5.6) 


5.5 


1688 


11251 


8.3 


4.9 


2.2E5 


7.0E2 


1.6E2 


41.3 


38.6 


M5 (10.0) 


5.0 


3000 


20000 


4.7 


4.4 


7.1E4 


2.2E2 


50.0 


40.8 


38.1 


3M4 (17.7) 


4.5 


5337 


35578 


2.6 


3.9 


2.2E4 


70.0 


15.0 


40.3 


37.6 


M4 (31.6) 


4.0 


9487 


63246 


1.5 


3.4 


7.1E3 


22.0 


5.0 


39.9 


37.2 


3M3 (56.2) 


3.5 


16876 


112509 


0.8 


2.9 


2.2E3 


7.0 


1.0 


39.4 


36.7 


Ensemble (127) 


7.27 








6.63 


1.57E7 


4.9E4 


1.1E4 


43.1 





Note. — All data listed in the above table are for single SC, except for the Ensemble. 
''GMC type (number of GMC in an ensemble). 
^GMC mass. 

'^Average gas density of a GMC. 
"^GMC core density. 
"GMC radius. 
'^Star cluster mass. 

^Star number for different stellar mass ranges. 

''Stellaj: bolometric luminosity. 

'Mechanical luminosity (Stellar wind + SN). 
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Bally et al. (1988) also obtained an average molecular gas density of 50 cm~^ for the region 
within a radius of 500 pc of the center of our Galaxy. In this paper, we assume a uniform 
ambient ISM with similar density surrounding the GMCs for our model. 

To investigate whether this medium should be considered as atomic or molecular, and 
to get an estimate of its mean density, we compare the observational constraints for various 
ambient ISM constituents (H2, HI, and H II) for the central 1 kpc region of M 82. Table [2] 
shows that gas with a column density of about 10^^ cm~^ is required by observations, and that 
the dominant state of the ISM is molecular. Hence, from this observed H2 column density, 
and an adopted diameter of 1 kpc for the starburst region, we can derive the number H2 
density of about 30 cm~^ (or a total H2 mass of ~ 1.2 x 10^ M©), and we adopt this figure 
for modeling the central 1 kpc region in M 82. 

In reality, intercloud gas i n M 82 is unlikely to be u niformly distributed, as assumed in 



our model. Recent studies (e.g. lGlover fc Mac Lowll2007l . and references therein) show that a 



smoothly distributed turbulent medium consisting of atomic gas would quickly (within a few 
10^ yr) develop density fluctuations, becoming a highly non-uniform medium of molecular 
H2, with density enhancements up to a factor of 100 or more times the mean density. For 
simplicity, we ignore these density fluctuations, and regard this medium as represented by 
its mean density, treating it as uniform for the purpose of computing the material swept-up 
by the expanding shells. 

The standard input parameters in our time- dependent PDR simulation that describe the 
physical properties of the model shell were summarized in Table 2 of Paper I. The cosmic-ray 
ionization rate ( adopted in this paper (the standard Galactic value, i.e. 1. 3 - 2.0 x 10~ ^^ 



s ^) is up to two orders of magnitude lower than that measured in M 82 (e.g. iFarquhar et al. 



I994J ). However, the main heating mechanisms of the PDRs are photoelectric, H2 FUV 
pumping, and H2 photodissociation. The contribution to the total heating rate from cosmic- 
ray ionization and excitation and the decay of turbulence within the cloud/shell is generally 

Table 2. Observed column density of various ISM components in the center of M82. 



Type 


Column Density (cm ^) 


Reference 


H2 


6.1 X 10^2 


Wild et al. (1992) 




a few X 10^^ 


Mao et al. (2000) 


HI 


2.6 X 10^2 


Weliachew et al. (1984) 


H II 


9.0 X 10^2 


Carlstrom & Kronberg (1991) 



negligible. These processes only becom e important at large depths during late starburst 
evolution and/or dark cloud (IBelll 120061 ) . 



The input parameters for a single-point dense dark-cloud model, which is used to pro- 
duce the initial abundances at the first time step (t = yr) all depth steps for our PDR 
simulation, were also described in Paper I. The dark-cloud assumption of chemistry is the 
same for all depth steps; it is a reasonable guess for the initial gas conditions in a GMC 
before star formation occurs. Our starburst model time step begins at t = 1.0 x 10^ yr, 
adopted as the time when the massive star formation occurs in the center of the GMC. For 
this time step and the subsequent time steps, the input abundances are re-set to the output 
abundances of the previous time step generated by the UCL_PDR code. The chemistry at 
first iteration is calculated from gas temperature and attenuated FUV flux, and then revised 
iteratively until the balance criteria of heating and cooling is reached for each depth step 
at each time step. Hence, the final results are not significantly dependent on our initial 
dark cloud chemistry input at t = yr. The metallicity dependence appears in several key 
processes in the our PDR calculations, and accordingly we adopt solar abundances for the 
metals, i.e. unit metallicity. The dust-to-gas mass ratio is adopted as 1/100. Standard values 
of dust properties are used in the model (see Table 2.3 in Bell 2006), though the UCL_PDR 
code allows the various dust properties to be specified as free parameters which can vary 
with shell depth and time. 

The chemistry within the parent GMC outside the shell is also handled by the same 
PDR analysis, using the different (lower) density in this region. The incident FUV strength 
for the cloud region is the attenuated radiation field emerging from the outer boundary of 
the shell, and the FUV strength inside the cloud is computed in the same way as for the 
shell, with the computation of Ay taking account of the lower density of the dust. 

Our multi-level non-LTE radiative transfer SMMOL code includes an empirical dust 
extinction model (see Table 1 in Mathis 1990). In our models, we assume the entire region 
containing all of the shells is unresolved. The lowest ten energy levels are incorporated for 
molecular species (CO, HCN, HCO+, CN, HNC), three levels for atomic [C I] and [O I], and 
two levels for atomic [C II]. In Paper I, single coUisional partner (H2) was used in the radiative 
transfer calculations. In this paper, multiple collisional partners (H, e~, 11+ , P-II2, 0-H2, He) 
are taken into account in the statistical equilibrium equation calculation. The collisional 
excitation of molecular lines involves two partners, i.e. P-H2 and 0-H2, but the excitation 
of [C I] fine structure lines is affected by collisions with all six particles, five (without He) 
for [O I] lines, and four (without He and H+) for [C II] lines. Since these forbidden lines 
have very low radiative transition probabilities, the upper states are populated primarily by 
collisions, and they are usually optically thin. 
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The input parameters to the SMMOL model are (1) molecular data including molecular 
mass, energy levels, transition frequencies, radiative rates and coUisional rates; and (2) 
physical data describing the object to model. This includes the physical distance of the 
current grid point to the center of the shell, gas density, number densities of the six collisional 
partners (H, e~, H"*", P-H2, 0-H2, He), the fractional abundance of molecules or atoms, the gas 
(kinetic) and dust (thermal) temperatures, shell expansion velocity, and the microturbulent 
velocity. 

Table [3] summarizes the parameters and variables used in our simulations. 

3. Model Results 

In this section we present examples of our simulation results for the shell ensemble 
centrally illuminated by massive star clusters. Detailed results for individual expanding 
shells can be found in Yao et al. (2006) and Yao Thesis. These shells are modeled in 
a similar way for both Winds and post-SN phases as described in previous sections. A 
family of these evolving shells form the basis of our starburst models in accordance with our 
description in § [2l Applications of the shell ensemble to M 82 and more distant starburst 
galaxies will be presented in subsequent sections. 

The two modeling phases are indicated by Winds and post-SN labels in tables and plots 
throughout the remainder of this paper. 

3.1. Kinematics of The Swept-up Gas 

The strong stellar winds and supernova explosions from hundreds to thousands of the 
massive stars fuel the hot bubbles over a timescale > 10 Myr. The kinetic energy in the 
supersonic wind is thermalized by a stand-off shock, and the high pressure downstream drives 
a strong shock into the ambient ISM. The swept-up gas condenses into a narrow shell as a 
result of radiative cooling. The wind mechanical luminosity Emech comes mainly from Wolf- 
Rayet (WR) stars, with some contribution from O stars. All other stars produce a negligible 
effect. 

During the Winds phase, the sizes of the initial Stromgren spheres in our model ensemble 
increase slowly with time. The Stromgren radius ranges from 0.02 to 4.9 pc with the number 
of Lyman continuum photons between 1.5 x 10^^ and 5 x 10^^ s~^ generated from the central 
star clusters derived from Equation (2) in Efstathiou et al. (2000). The wind bubble catches 
up with the ionization front of the compressed shell in a time less than 10^ yr. The strong 
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stellar winds cause the bubbles to expand quickly into their parent clouds and to sweep up 
more gas into the shells. When the most massive star in the most ma ssive star cluster (i.e. 



I2OM0 star in the M7 cloud) terminates as a supernova at ~ 0.8 Myr (IMac Low &: McCray 



19881 ). this marks the beginning of post-SN phase. At this time, the largest thin shell (M7) 
caused by the stellar winds is expanding at a speed of ~ 50 km s~^, and all the shells have 
swept up the material in their parent clouds. The Winds phase ends earlier (< 0.8 Myr) for 
shells smaller than that for the M7 cloud. After 0.8 Myr, the shells begin to expand into 
a less dense uniform ambient ISM (i.e. 30 cm"'^). The mechanical energy produced by the 
first supernova and the subsequent ones re-energizes the shell formed in the Winds phase. 

The hot bubbles begin to cool at ~ 0.8 Myr for the 3M3 shell and ~ 7.5 Myr for the 
M7 shell. At this time, the radius and velocity of the M7 shell are about 270 pc and 24 km 
s~^, respectively. After this time, the superbubbles start to lose their internal pressure, and 
the shell expansion velocity decreases rapidly. When the shell velocity approaches the sound 
speed of the ambient ISM, the shells should stall and become thicker and less dense. The 
latter effect is not included in our model, since the external pressure of the ISM is ignored. 
It is clear that the lifetime of the progenitor GMCs may be short, but the birth of massive 
star clusters and their impact on the surrounding ISM is profound. 

In addition, we calculate the total amount of gas Mmodei in out model GMC/shell 
ensemble following the evolution of the shells, as shown in Fig. [2l The discontinuity seen at 
1 Myr is caused by the phase change {Winds to post-SN), in which the parent GMC mass 
contained in the shell is no longer taken into account after the shell sweeps up all material 
in its parent GMC. This mass will be used as a template or reference value to be scaled to 
the total H2 gas mass in a observed region of M 82 using a analysis for our model line 
SEDs, under the assumption that the line flux in the measured region is proportional to the 
total molecular gas mass (see § IHfor details). 



3.2. Thermal Properties and Chemistry of the PDRs 

Over the 100 Myr of shell evolution the total mechanical wind power produced by 
individual cluster varies greatly, for example, from 10'^^ - 10^*° erg s~^ for the most massive 
cluster in our model ensemble. In this paper, the mechanical power profile is used only 
for obtaining average values over each phase in order to compute the shell dynamics. The 
kinetic energy of the shells is between 10^^ ergs and 10^^ ergs, depending on cluster mass. 
The FUV radiation incident on the inner surface of the shells [Ay = 0) is a function of time 
(Go oc Rg^; see in Paper I). The Go value is in units of the Habing field (1.6 x 10~^ ergs 
cm~^ s~^) throughout this paper. This value decreases from about 10^ - 10^ (depending on 
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Fig. 2. — Plot of tlie total H2 mass swept-up by the shells and the mass remained in parent 
GMCs ( Winds phase only) as a function of time. 
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cluster mass) at the onset of massive star formation (i.e. t = yr) to between 10^ and 10^ 
respectively at 5 Myr when most of the massive stars (M^, > 30 Mq) reach the end of their 
hfetime. At t = 100 Myr, the Go values drop by 4 - 5 orders of magnitude. 

3.2.1. Density and Temperature 

The shell density, thickness, and temperature are calculated for both Winds and post-SN 
phases. Fig. [3] shows an example of simulation results of an M7 shell. The density value 
varies from 10^ to 10^ cm~^ and the thickness is between 10~^ and ~ 10 pc over a 100 Myr 
period, depending on cluster mass. The plateaus seen at the beginning of the Winds phase 
are due to small changes in the expansion velocity and shell temperature. Before the shell 
sweeps up all of the material in its parent cloud (t < 0.8 Myr), the shell density declines with 
increasing shell radius and decreasing shell velocity, and the shell thickness increases with 
time. The dense phase of the shells (10^ - 10^ cm~'^) is very short lived (between 10^ - 10^ 
yr). After the first supernova occurs (i.e. post-SN phase), the bubble continues expanding 
adiabatically into a lower density ambient ISM until a time tc (indicated in the plots), when 
this hot interior begins to cool and the shell enters the snow-plow phase while conserving its 
total momentum. The shell velocity then decreases rapidly with a corresponding decrease 
in shell density, to about three orders of magnitude lower than that at the adiabatic phase. 
Such large variation in the shell density is due to the range of dynamic pressure produced 
by the range in the shell expansion speed. The shell thickness increases from 0.1 pc at the 
beginning of the post-SN phase to 10 pc at 100 Myr. Similarly, the thickness covers a large 
range because in the early phases, the shells are highly compressed and contain very little 
mass, so they are thin compared to later phases where these conditions are reversed. 

In the plot, the first big jump occurs when the wind shock front catches up with the 
ionization front, and the expansion changes from H II to wind driven. For smaller GMCs, 
this transition takes place in less than 10^ yr. The discontinuity (or gap) between Winds and 
post-SN phases is due to the model change from Winds to post-SN phase. A smaller jump 
is also seen when radiative cooling inside the bubble becomes dominant, the shell switches 
from the adiabatic to the zero pressure snow-plow phase. 

Fig. m shows an example of the gas and dust temperatures as a function of the visual 
extinction Ay for an M7 shell at different starburst ages. The Ay is set to be at the inner 
surface of the shell (i.e. boundary between the hot bubble and the shell), and increases 
toward the outer edge of the shell (i.e. boundary between the shell and its parent cloud 
or the ambient ISM). During the Winds phase, the cloud Ay progresses from the outer 
edge of the shell to the outer edge of the GMC (i.e. boundary between the GMC and its 
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Fig. 3. — Plot of the density [ug, solid line) and thickness {dg, dashed line) of an M7 shell 
as a function of time. The Winds phase is indicated by red curves, while the post-SN phase 
is indicated by blue curves. The radiative cooling of the hot interior occurs at tc indicated 
by the dotted lines. 
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ambient ISM). The gas temperature has a negative gradient from the inner edges of the shells 
to the outer edges, because the FUV flux is attenuated owing to dust extinction resulting 
in decreasing photoelectric heating across the shells. The FUV field strength G{t, Ay) at 
different Ay (or depth in the shell) is a factor of e~^-^^"^^ less than the flux at the surface 
of the PDR (or Gq). For example, at 1 Myr the FUV field strength at Ay = 2 (layer of 
C+/C/CO transition) is attenuated to ~ 6% of the value at the surface (Go ~ 10'^ - 10^) 
for the shells in the ensemble. The gas temperature is in the range 10 - 1000 K across the 
shells. It is about 1-2 orders of magnitude higher than the dust temperature at the surface 
of the PDRs. Fig. [5] shows an example of the temperature structure for an M7 cloud before 
the shell sweeps up all of its materials. The parent cloud is also heated by FUV radiation 
from the central star cluster. The minimum Ay for the GMC corresponds to the extinction 
due to the shell at the shell-cloud interface, and the maximum Ay is the extinction at the 
outer edge of the cloud. The gas temperature changes from 1000 K to about 10 K across 
the clouds, it drops rapidly with Ay once the FUV field is attenuated, reaches a minimum 
when the CO cooling is most efficient. But when the gas becomes optically thick toward the 
outer edge of the parent GMC, the gas temperature rises slightly with increasing Ay (e.g. 
age beyond 0.3 Myr) due to the assumption of a semi-infinite slab used in our PDR model. 
Since the cooling line emission does not escape from the outer edge of the GMC, the cooling 
of the gas is inhibited. This causes an increase in the equilibrium gas temperature at large 
Ay. Although this effect is nonphysical, it affects only our Winds model, and it has no effect 
on the post-SN phase, and thereby no effect on our derived age discussed later in § HI 

Although the physical properties of each giant molecular cloud and the star cluster 
born in its center vary greatly with cluster and cloud mass, the model profiles for the shell 
density, thickness and temperature are similar. This implies that different initial cloud 
conditions in a starburst environment may yield similar gas properties through the entire 
evolution. However, for individual shells, the physical properties of gas inside the shells 
change drastically with time. These gas properties that contain the imprint of different 
evolutionary phases, also determine the molecular line radiative transfer, and hence the 
spectral energy distribution of line fluxes. It allows us in principle to date the burst age by 
modeling the line spectrum energy distribution for various molecular tracers and comparing 
them with the observations of a starburst galaxy. 

3.2.2. Chemical Evolution 

The chemical structure inside the shell is stratified. The FUV photons are gradually 
absorbed and lead to relatively sharp transitions. In Fig. [6|, the transitions of atomic species 
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Fig. 4. — Plot of the time- dependent gas and dust temperatures as a function of visual 
extinction Ay for a M7 shell. Solid lines represent gas temperature, and dashed lines indicate 
dust temperature. The Winds phase model is indicated by red curves, and the post-SN phase 
model is indicated by blue curves. 
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Fig. 5. — Plot of the time-dependent gas and dust temperatures in the parent GMC (M7 
cloud) as a function of visual extinction Ay. Solid lines represent gas temperature, and 
dashed hues indicate dust temperature. 
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(H"^/H, C^/ C, O) to molecular gas (H2 and CO) in an M7 shell are shown. The H2 abundance 
becomes much more enhanced at Ay > 1, and the formation of CO occurs at Ay = 3-4. At 
the surfaces of the shells, the dominant coolant is the [O I] 63 fim fine-structure line. Deeper 
into the shells and the clouds the cooling by [C II] 158 /xm, [C I] 610 fim, and CO becomes 
dominant (not shown). The chemical structure inside the shells changes significantly for the 
first few million years. This further justifies the use of a time-dependent PDR model for our 
shell evolutionary models. Figs. [7] shows the chemical evolution inside an M7 cloud, before 
the shell sweeps up all of its material. 

3.3. FIR/sub- mm/mm Line Emission in Individual Shells and A Shell 

Ensemble 

3.3.1. CO and Its Isotope CO 

Fig. [8] shows our model ^^CO line SEDs (J = 1. . .9) as a function of the starburst age 
for a single shell (corresponding to the M7 CMC) and a shell ensemble. In plot (a) the 
total line flux Sco is summed from gas in one single expanding shell and its parent CMC 
with a mass of 10^ Mq (M7 CMC and Shell or SS model). The parent GMCs contribute 
significantly to the total lower-J line emission during the Winds phase. Our model results 
show during earlier Winds phase about 50% - 100% of total the ^^CO(l-O) line emission 
comes from the M7 cloud, but it decreases to less than 24% at the ^^CO(5-4) line, and to 
almost no contribution at J > 5. Our model M7 cloud has lower density than lower mass 
clouds (i.e. 3M3 - 3M6) and hence is less effective at exciting higher J transitions. The 
nonphysical small rise in Tgas at high Ay (see Fig. |5]) has almost no effect on the CO (1-0) 
line flux. At around 1 Myr, the line intensity drops by three orders of magnitude because we 
have not included the gas swept up in the GMCs in the subsequent model of the shells (i.e. in 
the post-SN phase). These model line fluxes predicted for individual expanding shells can be 
used as a comparison with future observations, for example, the known expanding supershell 
centered around SNR 41.9 -(- 58 in M 82, in order to constrain the physical conditions of the 
gas and the age of individual shells. 

In plot (b) of Fig. [S] the total line flux Sco is calculated from all shells and their parent 
clouds in an ensemble with 3.1 x 10^ < Mgmc < 10^ M© (Shell + CMC Ensemble or SGE 
model). If multiple transitions CO data for individual expanding shells become available in 
the near future, models presented in plot (a) could be useful to constrain the burst age and 
gas mass in the shell, such as the supershell described in previous section. More than 80% 
of the ^^CO line emission arises from the massive shells (3M5 - M7) in the ensemble. The 
line SEDs have two distinct maxima with one near the J = 6 - 5 transition and another near 
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Fig. 6. — Plots of the time- dependent chemical abundances of the main species (H, H2, 
H+, e~, C, C^, O, and CO) relative to the total hydrogen density, as a function of visual 
extinction Ay for an M7 shell. 
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Fig. 7. — Plot of the time-dependent chemical abundances of the main species (H, H2, 
H+, e~, C, C^, O, and CO) relative to the total hydrogen density, as a function of visual 
extinction Ay for the most massive GMC M7 in the ensemble. 
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the J = 3 - 2 transition. The first maximum is associated with burst age between 0.3 and 7 
Myr, and the second maximum is mainly associated with age older than 7 Myr. At age 0.2 
Myr, the two maxima (4-3,8-7) seen in the line SEDs are due to the sum of line emission 
of gas in the shells and parent clouds. It is clear that the CO excitation in the line SEDs 
varies with shell expansion or starburst ages. At ~ 1 (± 0.2) Myr ( Winds and post-SN phase 
transition), the Sco is a few orders of magnitude lower than those for other ages. This is 
an artifact of the switch from Winds to post-SN phase, where the GMC mass swept up in 
the Winds phase is not carried forward into the post-SN phase, and the continuity equation 
(or mass conservation) is applied to the less dense ISM (i.e. riism = 30 cm~^) instead of the 
GMC. Plot that shows the model line SEDs of ^^CO, HCN, HCO+, and Table that calculates 
the fraction of ^^CO, HCN, and HCO"'" line emission from individual shells and their parent 
clouds are presented in Yao Thesis. 



3.3.2. Atomic C, O, and 

The atomic forbidden transitions are the most important cooling lines arising in PDRs. 
The ratios of these lines and CO cooling lines can be used to derive the physical conditions 
in PDRs: for example, the incident FUV fiux Gq, gas density n and temperature T^as? as 
well as the ratio of Gq/ti. By comparing these model line ratios with observations, we can 
constrain the physical properties of atomic gas within a detected region. 

Upper left panels of Fig. [9] show the model line fiux (in Jy km s~^) for shell/GMC 
ensemble as a function of time for the [C I] 370/xm line at 809.3 GHz and [C I] 609/im line 
at 492.1 GHz for the most massive shell/GMC and the shell ensemble. The [C I] 370/im 
transition has an excitation temperature of 63 K and a critical density of 3 x 10^ cm^^. Both 
values are higher than the [C I] 609/im transition (24 K and 5 x 10^ cm~^). The atomic 
line fiuxes along with the molecular line fiuxes that we presented here are for the template 
model of the star clusters and molecular H2 clouds. The values for the actual masses for M 
82 will be derived from a fit of the fiuxes of this template model to the observed fiuxes. The 
[C I] 370/im line fiux emitted in the shells is generally higher than that in the [C I] 609/xm 
line, whereas in the cooler less dense parent clouds seen during the Winds phase the [C I] 
370/im to [C I] 609/im ratio is close to one (as seen in both plots (a) and (b)). The massive 
shells are the dominant source for the neutral carbon line emission in the post-SN phase. 
The discontinuity seen in the plots (near 1 Myr) is a result of switching phase from Winds 
to post-SN as explained previously in connection with molecular emission. Both line fiuxes 
increases with time in the post-SN phase. The integrated fiuxes of atomic lines predicted by 
our model generally increases with time, because more gas has been swept into the shells. 
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Fig. 8. — Plots of model ^^CO line SEDs for two different configurations. The Winds phase 
models are indicated by red dotted lines {t < 0.7 Myr), while the post-SN phase models are 
indicated by blue solid lines (1 < t < 8 Myr) and black dashed lines (8 < t < 100 Myr). 
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Upper right panels of Fig. [9] show the model line fluxes as a function of time for [O I] 
63 yum line at 4744.8 GHz and [O I] 145 fim line at 2060.1 GHz. The intensity increases 
with time for both lines, and then levels off after 8 Myr due to a sub-thermal excitation, the 
critical densities of these two lines are above 10^ cm~^. The [O I] 63 fim line flux is clearly 
stronger than the [O I] 145 /xm line throughout the entire starburst evolution. The [O I] 63 
/zm transition has a lower excitation temperature (228 K) than the [O I] 145 /im transition 
(326 K), and hence it is easier to be excited. 

Lower left panels of Fig. [9] show the model line flux of [C II] 158yum line at 1900 GHz 
as a function of time. The [C II] 158/im line has an excitation temperature of 92 K and a 
critical density of 3 x 10^ cm~^. The flux of this line increases with time. 



4. Understanding of Molecular Gas and Starburst Ages in M 82 

In this section, we apply our evolving starburst models by comparisons to an expanding 
molecular supershell centered around the supernova remnant SNR 41.9 + 58 in the starburst 
galaxy M 82, and to the multiwavelength data of the central 1 kpc regions of M 82, in order 
to arrive at some conclusions about the nature of these two regions. 



4.1. The Supershell Surrounding SNR 41.9 + 58 



Observations have detected an expanding supershell ce ntered around the bright SNR 



41.9 + 58 in both molecular line and radio continuum (e.g. IWeiss et al.l Il999l : IWills et al. 



19991 ). This supershell has a diameter of ~ 130 pc, an expansion velocity of ~ 45 km s~^, 
and a mass of ~ 8 x 10^ M c^. The kinetic eri ergy of the observed supershell is estimated 
to be about 1.6 x 10^'^ ergs (IWeiss et al.lll999l ). The kinematic evid ence for the supershel l 
appears most readily in the ^^C0(1 - 0) position- velocity (PV) plot (iNeininger et al.lll998l ) 
as a depression on the west side of M 82, bounded by a feature emerging toward lower 
velocities and possibly blended with emission associated with gas following orbits in the bar 
potential. Neininger et al. (1998) conclude that the depression seen in the ^^CO(l-O) PV 
plot coincides with peaks in emission of [Ne II] and radio recombination lines, providing 
evidence that the void is populated by ionized gas inside the supershell. Seaquist et al. 
(2006) show that their PV plot reveals no depression in ^^CO J = 6 - 5 but instead find a 
region filled with ^^CO J = 6 - 5 emission that is not evident in the underlying ^^CO J = 1 - 
map. Their line ratio PV map is consistent with the appearance of the channel maps, which 
show emission in the shell region extending over a very broad range in velocity. Seaquist 
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Fig. 9. — Plots of model atomic line fluxes (C, O, C^) as a function of time. Upper left 
plots (a) and (b): the solid curves are the [C I] 609 fim lines, and the dashed curves are the 
[C I] 370 fim lines. Upper right plots (a) and (b): the solid curves are the O I] 63 /xm lines, 
and the dashed curves are the [O I] 145 fim lines. Lower left plots (a) and (b): model [C II] 
158/im line flux as a function of time. The red color indicates Winds model, and the blue 
color indicated post-SN model. 
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et al. (2006) conclude that the location of this supershell contains CO with higher than 
average excitation, together with the ionized gas. The cavity created by the supershell is 
not associated with promi nent emission in higher density tracer such as HCN and HCO"'" in 



their low-excitation lines (iBrouillet fc Schilkdll993l : ISeaquist et al.lll998l ). This implies that 



the higher state of excitation may be due to higher kinetic temperature. Besides the known 
expanding supershell centered around SNR 41.9 + 58, there is evidence for other shells 
having sizes from several tens of parsecs to more than 1 kiloparsec, and kinetic energies 



between ~ 10^° and 10^^ ergs (e.g. 


Lo et al. 


1987 




Garcia-Burillo et al. 


2001; 


2002; 


Bartel &Bietenholz 


2005; 


Bavet et al. 


20081). 





The comparison of the kinetics of our single shell model with the observed supershell 
in M 82 is summarized in Table 3 of Paper I. Our model results and the observations agree 
remarkably well. In this paper, we investigate the state of excitation of the molecular gas in 
the supershell relative to that of the surrounding CO emitting gas in M 82, by comparing 
the predicted line ratios in the shell to those in the surrounding gas. For the surrounding gas 
we use line ratios computed for the bulk of the disk molecular gas based on our forthcoming 
analysis of fitting our model for a shell ensemble to the observed line ratios for the central 
1 kpc (see § I4.2.2p . Fig. (TD] shows this for the line ratios of ^^CO high J transitions to the 
(1-0) transition (i.e. Ico/Ico(i-o), Ico in units of K km s~^) for the model supershell. The 
jump in the ratios seen at J = 3 to 5 results from the addition of line emission of M7 shell 
to that of its parent cloud, where the CMC contributes 5 - 45% to the total line emission 
for J < 4, but less than 0.2% for J > 4. The plot shows clearly that our model for the 
supershell (red dashed curve) predicts that its line SED exhibits a higher level of excitation 
than the surrounding emission within M 82 (represented by the adjacent curve). Thus, one 
can expect some excess emission at high excitation transitions in the supershell after the 
underlying low excitation is subtracted out. Our model results are in qualitative agreement 
with the observational evidence for higher than average excitation emission in the supershell 



(e.g. iNeininger et al.lll998l : ISeaquist et al.ll2006l ). When higher quality and more extensive 
data on the excitation become available, our model predictions can be useful in interpreting 
the observations. 



4.2. FIR/Sub-mm/mm Line Emission in The Central Region 

Here we use the model components described in § [3] to produce a fit of our model 
line spectral energy distribution to the observations of molecular gas in the central 1 kpc 
region. The purpose is to determine whether it is possible to model the FIR/sub-mm/mm 
line emission in a massive star-forming galaxy, and whether there is a relation between the 
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Fig. 10. — Plot of model line ratios of ^^CO high J transitions to the (1-0) transition (i.e. 
^co / Ico{i-o)^ Ico in units of K km s~^) as a function of rotational quantum number J 
for an expanding supershell (M7) at age 0.7 Myr {Winds phase, red dotted curve). For 
comparison, a similar plot is shown of the observed (magenta circles with errorbars) and 
modeled SED (blue dotted curve) of the central ~ 1 kpc region of the disk of M 82 to 
represent the background disk emission with lower excitation. For details of the latter model 
fit, see § KT^. 
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molecular gas properties and the age of the starburst (i.e. finding the age indicator), and to 
assess the overall impact of the starburst on the fine scale structure and physical conditions 
of the ISM in M 82. 

We treat the entire central 1 kpc as an evolving starburst region, which can be modeled 
by following the evolution of an ensemble of expanding shells and clouds. Hence, different 
gas chemistry scenarios can be simultaneously at play in the center of this galaxy. However, 
our model does not attempt to reproduce or model the geometrical distribution of shells in 
an actual starburst system. In any event this distribution is unknown since the individual 
shells are not observed. The total line emission is assumed to be represented by the sum of 
the emission from all the shells in the model ensemble, which will then be used to compare 
with the observed data to estimate the stellar mass, the total H2 mass, and the age of the 
associated starburst in the measured region. 



4-2.1. Observational Data 



The central concentration (~ 1 kpc) of molecular gas in M 82, which fe eds the strong star 



formation activity, has been studied by naany authors since the 1980s (e.g. lYoung &: Scoville 



1984 : Wild et al. 
2000; Weiss et al 



992 



Gusten et al 



2001 



Ward et al 



19931: IWeiss et al.lll999l : iMao et al.ll2000l : iPetitpas fc Wilson 



2OO3I ). Interesting results arise from these studies. For 



example, the observed CO line SED and line ratios can be reproduced by emission from low 
(n(H2) ^ 10^ cm~^) and high (n(H2) ~ io3-5-4.5 > 4Q -^-^ excit ation gas com- 

ponents using a Large Velocity Gradient (LVG) method (IWeiss et al.ll2005l . and references 
therein). The high excitation component, responsible for the excitation of levels beyond J 
= 4, arise from dense and warm gas, while the low excitation component is emitted by dif- 
fuse low density gas. The LVG method assumes a uniform abundance and velocity gradient 
across the modeling region, and no star formation history is considered as a cause for these 
conditions. It is this singular distinction which is the focus of this paper. 

The excitation conditions of multiple transitions of dense gas tracers HON and HCO^ 
in M 82 have also been investigated, for example, by Seaquist & Frayer (2000). It was found 
from an LVG model, that both species are excited under a common set of conditions in star- 
forming regions where the n(H2) is near 10^ cm~^, Tfcj„ = 50 K, and the abunda nces of HCN 
and HCO+ are 2 x 10~* and 1 x 10"^, respectively (jSeaquist fc Fraverll2000l ). Molecular 
lines are commonly observed at 22 beam size, which covers about 680 pc of the center with 
a total H2 mass of a few times 10^ M0 in M 82. 



The atomic coolant, far-infrared lines in M 82, e.g. [C I] 370 fj,m, 609 /zm, [O I] 63 fim, 
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146 Atm, and 



Colbert et al. 



C II] 158 Atm, have been stud i ed by several groups (e.g. IStutzki et al.lll997 



19991 : IPetitpas fc WilsonI l200ll : iLord et al.lll996l : iNerishi et al. 



2001 



and ref- 
erences therein). As is the case for the molecular lines, the ratios of these cooling lines may 
be used to constrain physical parameters and possibly the age of the starburst. These ratios 
are sensitive to the physical and chemical conditions (density, temperature, and abundance), 
hence provide an opportunity to model the physical state of the neutral gas. In addition, 
unlike optical atomic line tracers, these FIR lines are relatively insensitive to extinction. A 
close examination of these lines emitted in M 82 may provide a template for future com- 
parisons to infrared-bright, dust obscured starburst galaxies like M 82, including those at 
high-2. 

In order to provide a useful indication on the starburst age(s), it is desirable to make 
comparisons with multiple transitions for various molecules and atoms. However, meaningful 
comparisons can be made only for regions where observations refer to the same beam size. 
The diagnostic tracers used in this paper are molecular ^^CO, its isotope ^^CO, HCN, IICO+, 
and atomic C, O, and C^. The low- J ^^CO lines are easily excited at relatively low densities 
and temperature, and are found essentially in every molecular gas cloud, and so they are good 
diagnostic tools for total molecular H2 content, diffuse gas conditions, and star formation 
history. The less abundant ^^CO isotope has a much lower optical depth, and the line ratios 
between optically thin transitions in ^^CO are more reliable probes of the total gas content 
than ^^CO. The CO molecule is not considered a good tracer of dense and highly excited 
gas that is directly involved in starburst (i.e. earlier phase of star formation). However, 
molecular HCN and HCO"^ lines are more sensitive to dense gas (i.e. pre- or post-birth 
of stars) owing to their higher critical densities than CO. The atomic C, O, and C^ fine 
structure lines are excellent probes of the PDRs in starburst regions, and their line ratios 
can be used for diagnosing the conditions of the associated FUV flux and gas density, as well 
as for indicating the ages of the later stages of starbursts. 

In this paper, we use the observations of molecular and atomic gas in the central 1 kpc 
of M 82 described above. The molecular data for the center 1 kpc (along the major axis of 
the disk) are taken from Weiss et al. (2005). T he molecular data fo r the center 680 pc are 
taken from Mao et al. (2000). The atomic data (iNegishi et al.ll200ll . and references therein) 
to be used in our ratio-ratio diagram analysis, which are obtained from a larger area (~ 
1.2 kpc) than for the molecular data. Note that the atomic C data are not included in our 
l ater ratio-ratio ana lysis, because the two [C I] line data correspond to different beam sizes 
( IStutzki et al.lll997l ). These molecular data will be used in comparisons with our models. 
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4-2.2. Model Fit to the ^^CO Line Spectral Energy Distribution 

We consider first a model involving a single instantaneous starburst, and later consider 
whether extended starbursts could also provide an adequate fit. 

Part I: Instantaneous Starburst Model 

Our initial starburst model involves a single event in which all of the stars/clusters 
are formed simultaneously and instantaneously in the center of GMCs, associated with a 
unique age t and a star cluster mass (corresponding to a GMC mass Mgmc = 4M^, at 
the birth time). An instantaneous model, though physically unrealistic, is an acceptable 
representation of the SED if the duration of the star forming event is short compared to the 
age of the starburst. The intent is to derive these two parameters by fitting to the data. By 
extension, the total mass of H2 swept-up in the ISM at any age is also determined. 

The method used for fitting is the well-known numerical chi-squared (x^) procedure. 
The number of data points is six for the ^^CO data, and the number of free parameters 
(in this case) is two. The quantity S'^^^^i represents the corresponding model to be fitted, 
SLdeiifi t) = f Si^p (t), where Si^p{t) is the model template hue SED at age t, as given 
in Fig. [H] (see § [3]), corresponding to a model template GMC mass M^^^ , a model template 
cluster mass M^^^^, and a model template swept-up mass by the shells M/^'^^. The values 
for these parameters are Mg^^ = 1.69 x 10^ Mq, and M*^^^ = 4.2 x 10^ Mq. These initial 
masses correspond to the 127 clusters included in Tabled] (see §|2]). The ratio of stellar cluster 
to GMC mass is 0.25 according to the assumed SEE. The adjustable dimensionless parameter 
/ is introduced to control the amplitude of the model line SED (and hence the total cluster 
mass), and the age parameter t controls its shape and slope. These are simultaneously 
adjusted to provide the best fit corresponding to the minimum xl- By assumption, the line 
fluxes S^^^pit) are summed over the contributions of all clusters and GMCs, so that the best 
fit GMC mass Mgmc, cluster mass M*, and the shell swept-up mass Mg^ are determined 
from the corresponding best fit value of the parameter / by the relations, 



Mgmc — fM^^^^ , (3) 
M, = /M;_^, (4) 
MsK = fM,f^p (5) 

We calculate xl for a range of t and /. A minimum xl value is obtained with a standard 
error estimation (i.e. the traditional likelihood method from using an inverse Hessian matrix 
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or covariance matrix). In the results there were found to be two minima with acceptable 
values of x^(min), one for the Winds phase and the other for the post-SN phase. The chi- 
squared contour surrounding the minimum for the post-SN phase is shown in Fig. [TTl where 
the contours are xl = xli'"^^^) + ^ / - P)? and i = 1, 2, 3, . . ., corresponding to contour 
intervals of la. The existence of a numerically acceptable solution for each phase signifies 
that there exist conditions in the molecular clouds of the Winds phase which are similar 
to those found in the compressed shells associated with the post-SN phase. However, only 
the post-SN phase solution is acceptable physically, since an age of 0.07 Myr derived for 
the Winds phase is implausibly small for a variety of reasons. For example, it is impossibly 
short compared to the dynamical time for the region (a few Myr) which would control the 
duration of the starburst. Such a small age might be barely plausible for an individual shell, 
but not for the molecular gas occupying this entire region. 

The best fit line SEDs and ^^CO line SEDs at 4 and 7 Myr are shown superposed on the 
data in Fig. [121 The initial stellar mass and GMC mass are relatively small, but the impact 
on the surrounding ISM is significant. The total H2 mass swept up by the shells is ~ 2.0 ± 
0.1 X 10^ M0 at the best fit age of 5.6 Myr. This predicted value is in good agreement with 
i-u^ ^^r.r. +1 — „ — 1 region obtained by other studies (a few 10^ Mq) 



Rieke et al. 


1980; 


Mao et al. 


2000) 



(e.g. 

is ~ 1.4 X 10^ Lq using information from Starburst99 based on our instantaneous starburst 
model. 

In order to investigate how sensitive the results are to the assumed initial upper mass 
limit of the cluster spectrum (and corresponding GMC mass spectrum), we repeated the 
above analysis with revised upper mass limits of both 7.5 x 10^ M© and 2.5 x 10^ M© for 
the stellar spectrum, and corresponding GMC upper mass limits of 3 x 10^ Mq and 10^ Mq 
respectively. For the first case, we find t = 5.0 ± 0.4 Myr, / = 1.8 ± 0.2 with xl = 0.9, and 
for the second case an unacceptable fit with a. xl = 33. Thus, a comparable solution may 
be found with a choice of a slightly lower upper mass cutoff, but no acceptable solutions are 
found with values reduced by a factor of 10 or more in the upper cutoff of the cluster mass 
spectrum. We conclude that the model can provide acceptable fits to the data only if the 
dominant initiating starburst clusters are massive, at least 5 x 10^ Mq. 

Our models show that the H2 density of the shells at the best fitted age 5.6 Myr is 
between 10^ and 10^ cm~^, and the gas temperature is ~ 50 - 100 K . These values are 



comparable with the two-component LVG predictions (jWeiss et al.ll2005l ). The evolution of 
CO abundance as a function of Ay is illustrated in Fig. [HI The CO abundances in massive 
shells (M6 - M7) are above 10^^ with respect to the total H density, providing most of the 
CO emission. 




Fig. 11. — A contour plot of values as a function of mass coefficient / and burst age t at 
the post-SN phase. 



- 36 - 



M 82 (~1 kpc) 




I I I I I I I I I I I 

01 23456789 10 

'^CO Rotational Quantum Number J 



Fig. 12. — A fit of an instantaneous starburst model to the ^^CO line SED for the central 
1 kpc disk region of M 82. The blue curve is the best age at 5.6 Myr for the post-SN phase, 
the cyan curve is the line SED at 4 Myr, the green curve is the line S EP at 7 Myr, and the 
observed data are indicated by magenta open circles with error bars (IWeiss et al.l 120051 ). 
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Part II: Extended Starburst Model 

The foregoing discussion and results assume an instantaneous starburst with the result 
that our CO best fit model has an age of 5.6 Myr. The question naturally arises whether 
a model with a period of more continuous star formation would also provide a satisfactory 
solution. One can anticipate that the answer might be yes, if the best fit line SED were to 
be roughly equally represented by the SED of an outburst at one epoch, or an average SED 
over some time period roughly centered on, and symmetrically distributed about this epoch. 
Fortunately, it is straight forward to test this hypothesis since the SED for a smoothly varying 
star formation model may be constructed from a superposition of instantaneous bursts at 
different times. 

We are thus led to consider the extreme case of a uniform star formation rate (or 
SFR) occurring between an epoch 10 Myr ago and the present time. This starting point of 
the event may be considered appropriate because the most massive shells from even earlier 
epochs would now be large enough to exceed the thickness of the nuclear disk and thus their 
emission would begin to fall outside the region modeled. Fortuitously, this period is also 
almost symmetrically distributed about the epoch for the best fit instantaneous model. 

We do not discuss the procedure in detail here, since the analysis proceeds as before, 
but with only one parameter, namely the star formation rate over the past 10 Myr. The 
SED employed is then an integral of the CO line SED profiles over a time period of 10 Myr. 
The result is that an acceptable fit (minimum xl = 0.78) can be found for a continuous SFR 
— 0.5 ± 0.05 M0 yr~^. The total stellar mass produced during this period is (5.0 ± 0.5) x 
10^ Mq, which is, not surprisingly, close to the total mass (i.e. 4.3 x 10^ Mq) required for 
the single epoch model. The conclusion is that a uniform SFR over the past 10 Myr also 
produces a satisfactory fit to the -"^^CO data. In addition, it may be plausibly inferred that 
a variety of star formation histories would work, provided the SFR rate profile is more or 
less symmetrically distributed about the epoch of 5.6 Myr. 

The implication of this result is that in terms of the agreement between the model 
and the data, the star formation need not be instantaneous, or even sharply peaked at 5.6 
Myr. However this epoch nevertheless represents a unique point of time associated with 
the history of star formation in M 82 since it would emerge from various representations of 
the star formation profile. If there were an instantaneous starburst in M 82, then its age 
is approximately 5-6 Myr (plus or minus a factor of about 2). However, the SED is also 
consistent with continuous star formation rate of 0.5 Mq yr~^ over at least the last 10 Myr. 
The factor of 2 is our estimate of the uncertainties introduced by the uncertain density of 
the ambient interstellar medium in the model, as well as our assumption that the shells do 
not stall. 
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4-2.3. Molecular and Atomic Line Ratio Diagrams 

The line intensity/flux ratio-ratio diagram can be another diagnostic tool for studying 
the gas excitation conditions and properties, as well as their relations to starburst evolution, 
especially when there are not enough data points available for the type of model fitting 
discussed in the previous section. Since the line ratio is independent of the total gas mass in 
the measured region, the ratio-ratio diagram cannot be used to provide an estimate of the 
total gas mass directly. However, once the age t is obtained, we can scale the template line 
flux spectrum to the flux observed, and calculate the model value for the swept-up H2 mass. 

12 CO and Its Isotope ^^CO 

Fig. [13] illustrates the ratio-ratio diagrams for different transitions involving ^^CO and 
^^CO predicted by our model {post-SN), and a comparison with the observations. The 
observed data refer to the center of M 82 with a beam-width of 22 (see Table 1 in Mao et 
al. 2000). All line brightnesses are compared in units of Jy km s~^. The isotope abundance 
ratio [^^CO]/[^^CO] of 55 is adopted for the ensemble modeling. In the plots, we include 
the systematic uncertainties (31% for ^^CO{7-6), 20% for ^^CO{4-3), 16% for ^200(3-2), 
23% for ^^00(2-1), 10% for ^^00(1-0)) into the line ratio error estimations (i.e. sizes of 
error bars). In plot (a) the model ratios of ^2CO(7-6)/(4-3) versus (2-l)/(l-0) match nicely 
(as expected) with the observations at age 5-6 Myr for the center 680 pc region. It is 
similar to t he age derived fro m the chi-squared fit to the ^^CO line SED in the center 1 



kpc region (jWeiss et al.ll2005l ) even though the angular size of the region is different. This 
is expected, since ~ 65 - 80% of the ^^CO emission from the inner 1 kpc disk originates 
from the central 680 pc starburst regions. Plots (b) and (c) show a poor match between 
our model line ratios of i3CO(3-2)/(2-l) versus (2-l)/(l-0) and i2cO(2-1)/^3qo(2-1) versus 
^2CO(1-0)/^^CO(1-0) (blue-dashed curves) and the observed data. The closest match within 
the observed uncertainties is 7 - 8 Myr for plot (b), 5 - 6 Myr for plot (c), where the latter 
is in a fair agreement with ^^CO best fitted age. 

Our model fails to produce the right ratios for lines involving ^^CO. If the choice of 
the isotope abundance ratio is to be considered as the reason for such poor fit, adopting a 
different isotope abundance ratio (55 is used in this paper) can affect the result in plot (c) 
but not that in plot (b). Mao et al. (2000) indicated that their ^^C0(2-l) values should be 
considered with caution, due to uncertainty of convolving a smaller beam (13 ) to a larger 
beam size (22 ). If we assume that the best match age for plot (b) and (c) should be between 
5 and 6 Myr, and if we assume that an erroneous value for the ^^C0(2-l) model flux is the 
reason for lower ratios seen in plot (b) and (c), we estimate that this value is underestimated 
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Fig. 13. — The ratio-ratio diagrams of molecular ^^CO and ^^CO line intensities (in units of 
Jy km s~^). Model results for a shell ensemble are indicated by the crosses connected with 
blue dashed lines. The age sequence is 2, 3, 4, 5, 6, 7, 8, 9, 10, 20, 30, 40, 60, and 80 Myr. 
The magenta filled circles with error bars are the observed data (22 resolution data for the 
center of M 82 from Table 1 of Mao et al. 2000; the errors include systematic uncertainties). 
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by a factor of 1.5. Another factor that contributes to the poor match between our model 
results and the data is that the line ratios involving ^'^CO may be particularly sensitive to 
optical depth if the intensities are not optically thick. 

We also investigated the effect of reducing the upper mass limit to the GMC mass 
spectrum as was done for the ^^CO model, and found that changing the upper mass limit to 
the GMC mass spectrum has no effect in resolving this problem. 



Atomic and 



Fig. [H] shows the model ratio-ratio diagram for [O I]63/im/[C II]158//m versus [O 
I]63/im/[0 I]145/im, and a compari son with the observ ations of these atomic lines from 



the central 1.2 kpc region in M 82 (INegishi et al.ll200ll ). All line fluxes are compared in 



units of W mr"^. The model [C II]158/im line flux may be underestimated, since we ignore 

the line emission that arises from the H II region. A good match between our model and 

the observation is obtained with age t ~ 10 Myr. The age predicted from atomic data is 

older than the age (~ 5 - 6 Myr) derived from our ^^CO line SED analysis. This may be 

// 

because the atomic line data are based on a 80 x 80" beam area whereas the CO line data 
pertain only to the 60 x 18 beam area. We suggest that these two ages may be a result 
of sampling different regions. More discussion of this possibility will be given in § 14.3.21 

The ranges of gas conditions for the model shells at 10 Myr are Gq ~ 350 - 1.4 x 10^, 
n(H2) ~ 10^ - 2.4 X 10^ cm~^, and Tgas > 20 K. The gas conditions derived from our atomic 
models for this sampling region are comparable with the study by Colbert et. al. (1999) (Go 
= 630, n = 2.0 X 10'^ cm~^), but the age is greater than that (3-5 Myr) derived by Colbert 
et al.. 

The total molecular gas swept up into shells cannot be obtained directly from the ratio- 
ratio diagram. However we can obtain this from the ratio of observed [O I] 63 /im line flux 
to the model template [O I]63 /im flux at age 10 Myr, i.e. / = Sobs / Sl^^J^ei ~ ^-^^ > where 
Sobs = 169 X 10-^5 W m-2 and S^^^^i = 202 x 10"^^ W m'^. We compute the total H2 gas 
in the measured 80 region by multiplying / = 0.84 by the model template H2 mass M^^^^. 
Hence, the resuh is M(H2) at age 10 Myr is ~ 3.4 x 10*^ Mq. 
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Fig. 14. — The ratio-ratio diagram of atomic fine structure line fluxes (in units of W m ^). 
The models are indicated by open circles connected with a blue-dashed curve for post-SN 
phase. The age sequence is 1, 2, 3, 4, 5, 7, 8, 9, 10, 20, 30, 40, 60, and 80 Myr for the post-SN 
pha se. The filled circle with error bars show the observed data for the center 1.2 kpc of M 
82 jNedshi et al.ll200lh . 
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4.3. Discussion 

4-3.1. An Expanding Supershell Associated with SNR 41-9+58 

The very good agreement between our supershell kinematic model and the observations 
is consistent with the hypothesis that this expanding supershell is created by strong mechan- 
ical winds from a young star cluster with a total mass of about 2.5 x 10^ Mq which formed 
at its center about 0.8 Myr ago. This agreement also suggests that the set of models we have 
put forward in this paper may be used to interpret other shells in M 82 or shells in other 
starburst galaxies. Although hke any other model, the result depends to some degree on the 
set of initial conditions and assumptions that we selected for our models. The rehability of 
the age and mass for this supershell derived from our kinematic study needs to be further 
examined in near future when high resolution maps of multiple transitions of CO emission in 
this shell are available to compare with our model. Meanwhile we relate our model CO line 
ratio SED (i.e. Ico/Ico{i-o) as a function of J) for the SNR 41.9 + 58 (i.e. M7) supershell 
at age 0.8 Myr to the corresponding line SED for the surrounding gas in the inner 1 kpc 
starburst region of M 82. We show that the emission in the M7 supershell exhibits a higher 
degree of excitation than the surrounding emission (see Fig. [TOl) . 

There are a number of issues arising from the supershell study. They are as follows: (1) 
it is interesting to ask whether our results are consistent with a possible physical association 
between the supershell and the bright SNR 41.9 + 58 near its center. If the bright SNR were 
within or near the SSC, there may not be sufficient gas remaining to form an SNR after the 
action of the winds from the cluster; and (2) the SSC responsible for the formation of the 
supershell might also have provided the stellar mass for the several h undred solar mass blac k 



hole detected by Chandra X-ray observations near its center (e.g. iDewangan et al.ll2006l ). 
Theories for the formation of this black hole include the collapse of a hyperstar formed by 
the coalescence of many norm al stars, or the direct merger of stellar mass black holes (e.g. 



Kawakatu fc Umemural l2005l ). The SSC is adequately endowed with sufficient mass since 
there would hav e been 1,700 O stars, each with mass greater than or equal to about 40 Mq 



(lYao et al.ll2006h . 



4.3.2. Central Starburst Region 
Age of Recent Starburst and Star Formation History 



Given the complexity of M 82, a full understanding of star formation epochs requires 
various diagnostic tools to trace different ISM components in starburst regions. Especially 
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since both optical and near- to mid-infrared emissions are subject to higher extinction in 
dusty media, the selection of SSCs may be biased toward either younger or older age as 
mentioned before. Since the ISM is nearly transparent to FIR/sub-mm/mm emission, the 
analysis in this paper, employing atoms and molecules emitting in this range, forms a useful 
complementary investigation to those already mentioned. Using our evolving starburst model 
for neutral gas media, we have been able to probe the recent star formation history of M 82 
throughout the entire volume of the central starburst region. 

The starburst ages derived from optical and infrared spectra are 5-6 Myr, 10 - 25 
Myr, and 30 - 100 Myr. The age derived from our analysis of CO line SEDs and ratio-ratio 
diagrams is 5 - 6 Myr for the central 1 kpc x 280 pc rectangular regions; although the region 
used is 1 kpc, about 70% is concentrated toward more central regions (~ 350 pc). The age 
derived from our atomic data is slightly older, i.e. 10 Myr for a larger area (~ 1.2 kpc). 
We suggest that these two ages may be a result of sampling different regions as mentioned 
earlier. It is unclear from our analysis whether these two ages refer to the same period of 
star forming activity or to two spatially separated independent bursts. A more sophisticated 
model and more data are needed to clarify the picture. 

The burst ages derived from our model are similar to the results found in the aforemen- 
tioned studies by Forster-Schreiber et al. (2003) and Efstathiou et al. (2000). However, for 
the atomic data there is a discrepancy between our result (10 Myr) and the study by Colbert 
et al. (1999) (3 - 5 Myr) using a similar set of data (by Negishi et al. 2001). Nevertheless, 
our derived gas conditions for the shells at 10 Myr {Gq ~ 350 - 1.4 x 10^, n(H2) ~ 10^ - 
2.4 X 10^ cm~^, and Tgas > 20 K) are similar to those derived by Colbert et al.. The age 
discrepancy may be caused by differences in the choice of models, for example, Colbert et 
al. used the CLOUDY PDR model to compute the atomic line fluxes, while we use the 
UCL_PDR model and a non-LTE radiative transfer model to compute the line fluxes. 



Molecular Gas Properties 



Our evolving shell models yield familiar values for the gas density, temperature, and 
structure scales compared to those measured in the center of M 82 fe.g. iLynds fc Sandage _ 



1963 



2001 



Rieu et al.l 



Ward et al. 



989l:IStutzki et al.lll997l : ISeaquist fc Fravei]l2000l : iMao et al.ll2000l : iNegishi et al. 



20031 ). For an extended starburst scenario {SFR = 0.5 ± 0.05 M© yr~^ as 
discussed in § I4.2.2p . the shell densities are in the range 10^ - 10^ cm~'^, and the gas temper- 
atures are in the range 20 K to 1000 K across the shell for various shells. The total H2 mass 
swept up by the shells within the inner 1 kpc (~ 2.0 ± 0.1 x 10^ M©) and 1.2 kpc (~ 3.4 ± 
0.3 X 10*^ Mq) detection regions are c ompatible with those derived from the CO luminosit y 
using the CO-to-Hs conversion factor Jwild et al.lll992l : iMao e-TaPboOol : Iwalter et al.ll2002h . 
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It is also comparable with the total ambient gas mass in our model. Hence, the picture sug- 
gested is that of a porous neutral ISM in the central star-forming region of M 82, a product 
of evolving shells. In reality, many or most shells are probably in the form of fragments, 
small cloud clumps, sheets, or partial and full circular arcs (e.g. 
20061 . and references therein). 



Lo et al. 1987: Yao et al. 



Molecular Ring Formation Mechanism 



Although different stages of starburst evolution are applicable to different central regions 
of M 82, the shell sizes and the physical conditions of the gas within the rings (diameter ~ 
300 - 600 pc) predicted by our model are similar to what is expected from models involving 
expanding shells from a central starburst such as those proposed by Carlstrom & Kronberg 
(1991). Their hypothesis is that molecular rings in M 82 are a result of compressed gas in a 
starburst region. This hypothesis is supported by the observations of the geometrical struc- 
ture of the CO line emission and continuum emission, as well as the discovery of supershells 
that have not yet had time to break out of the galactic plane. However, the conclusion drawn 
from the shell size and average gas conditions in the inner 1 kpc region is only suggestive, 
since our model does not handle the physical distribution of molecular gas in the center of 
M 82. It is also important to realize that the foregoing interpretation of the lobes as a ring 
is not unique. A number of authors have argued that the molecular rings are a product 

related with the gravitational effects of the bar (e.g. 



Shen & Lo 


1996; 


Wills et al. 


2000) 



4-3.3. Limitations of Our Model and Their Impacts 

We have demonstrated that the kinematic and FIR/sub-mm/mm emission properties 
of individual expanding shells and star-forming regions in a starburst galaxy like M 82 can 
be understood by following the evolution of individual massive super star clusters or an 
ensemble of such clusters surrounded by compressed shells and GMCs. It is an important 
piece of complementary work to the existing optical and infrared studies, and it helps us to 
obtain a more complete and or accurate picture of star formation episodes in the center of 
M 82. 

However, our model also has a number of caveats, limitations, and potential sources of 
systematic error. Here is the important list: 

(1) We have neglected throughout the effects of the ambient pressure in slowing down 
and perhaps stalling the shells. This applies to both Winds and post-SN phases. We recall 
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that the shells will stall when their expansion velocities decrease sufficiently that they are 
approximately equivalent to the sound speed (P/p)^ of the external medium. To estimate 
the effects of this pressure, we can thus compute the sound speed associated with estimates 
of the pressure and compare this with the shell speeds. We compute the total pressure 
Prinnd luslde the cloud, a ssuming it is in virial equilibrium from the following equations 
(iMcCrav fc Kafatoslll987l ). 



Pcloud Pext ~l~ Pinterni (6) 

Pexi = 2noA;T, (7) 

Pintern 

0.5GS' (8) 



where Pexter is the external pressure, Pmtem is the internal pressure, no is the cloud H2 
density, k is the Boltzmann constant (1.38 x 10^^^ erg K~^), G is the gravitational constant 
(6.67 X 10~^ cm^ g"^ s~^), and S = Mgmc I {t^ P-gmc)- The sound speed in a GMC can 
be calculated from the equation. 



Ccloud 



/ Pdoud \ 2 
^ Pdoud ^ 



(9) 



We obtain a sound speed of Cdoud = 19 km s~^ in a M7 cloud (n(H2) = 300 cm~^, Mgmc 
= 10^ Mm, and Rg mc = 47 pc), assuming Pext/k = 10^ K cm~^ in starburst regions of M 82 
(ISilich et al.l 120071 ). where k is the Boltzmann constant. We can further combine Equations 
(16) through (19) with the cloud relations Equations (7) and (8) to furthermore yield the 
sound speed in any given GMC, 



Cdoud = 19 kms ^ (^^^) * (10) 

For the external ambient medium, we use the aforementioned external pressure to obtain 
the sound speed for the ISM (n^sm = 30 cm~^). 



Cism = 40 kms (11) 

The comparison between sound speeds inside the clouds (Equation (20)) and shell ex- 
pansion velocities for the Winds phase indicates that shells from cloud masses above 10^ Mq 
would not be trapped, and those equal or below this mass would be stalled if the effects of 
cloud pressure were included. 
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We also compare the sound speeds in the ISM (40 km s~^) with shell expansion velocities 
for the post-SN phase. For example, for the shells associated with the three most massive 
GMCs, namely M7, 3M6, and M6 in our model ensemble, the shells have radii of 220, 180, 
and 130 pc at the best fit age 5.6 Myr without ambient ISM pressure, respectively. But the 
stall radii and ages for these three shells are 112 pc at 2 Myr, 70 pc at 1.2 Myr, and 42 pc 
at 0.8 Myr, respectively. Thus, without the inclusion of the effects of pressure it may be 
said that the shell radii at the time of observation are probably overestimated by more than 
a factor of two compared to the stall values when pressure is included. Since the swept-up 
mass by the shell is proportional to the R^, an overestimate by a factor of two in shell radius 
would yield a factor of eight in the total swept-up mass for a given GMC/SC mass. This 
may help to understand the shortfall in IR luminosity predicted by our starburst model (see 
point (3) for detailed discussion). 

Another issue worth mentioning is that the confining pressure will vary greatly with 
location in the galaxy, especially between the center and the edges of the disk where some of 
the observed supershells are located. For example, as we mentioned earlier in this chapter. 



Weiss et al. 


1999: 


Wills et al. 


1999) 



This supershell has a diameter of ~ 130 pc, an expansion velocity of ~ 45 km s~^, and a 
mass of ~ 8 X 10^ M0. If Pext/k = 10^ K cm~^ were the relevant external pressure in this 
case, then this shell will stall soon. However, the pressure may well be lower than the above 
value in this region, since part of the shell is seen outside the disk. Other expanding shells 
(incomplete arclike shapes) with velocities po ssibly as low a s 10 to 15 km s~^ with radii ~ 



200 pc are also observed in the central region (ILo et al.lll987l ). suggesting a sound speed less 
than the 40 km s~^ figure used above. 

(2) Observations of nearby bubbles in our own G alaxy and in the Magellanic Clouds indi- 
cate that the simple adiabatic bubble/shell theory (IWeaver er al.l 119771 : iMcCray fc Kafatos 
19871 ) coupled with the mechanical luminosities calculated by Starhurst99 for this paper 
leads to significant overest imates of the bubble pressure and hence the shell radius (e.g. 
Oey fc Garcia-Segurall2004l . and references therein). Either the wind power is lower or some 
hot bubble gas escapes from the bubble interior. In addition, Dopita et al. (2005) argued 
that the conventional bubble/shell dynamical model may overestimate the winds and super- 
nova mechanical power. Another argument is that gravitational instability may induce new 
star formation inside the shells. If such effects were present, they would have an impact 
on the estimate of the total stellar mass and luminosity in our model, as described in more 
detail in point (3). 



(3) The bolometric luminosity for the best fit cluster mass of 3.7 x 10® M© and best fit 
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age of 5.6 Myr is 1.5 x 10^ L0 (based on Starhurst99 model). The observed IR luminosity of 
M 82 disk is about 3.0 x 10^'^ Lq. Since the bolometric luminosity should be an upper bound 
to the IR luminosity from the same stars, the shortfall in the model luminosity is at least 
a factor of 20. This shortfall is similar to that (also about a factor of 20) between the star 
formation rate derived from our continuous star formation model (~ 0.5 yr~^), and the 



measured star formation rate ^ 5-10 yr ^ for the center of M 82 (e.g. Ide Grijs et al. 



2OOII : iLipscy fc PI avchani I2OO4I ). Hence, the stellar cluster mass needed according to the 
model to produce the observed CO luminosity is not sufficient to account for all of the 
stellar luminosity or young stellar mass in M 82. 

There are several reasons that our shell ensemble model may have overestimated the 
swept-up gas mass and the line emission for a given cluster mass, or equivalently underes- 
timated the stellar mass and luminosity required for a given swept-up gas mass. Points (1) 
and (2) above show that our model itself may be fundamentally optimistic in its impact on 
the ISM, i.e. the model shells may be too big for the stellar mass which generates them, thus 
leading to an overestimate in the swept-up shell masses and CO luminosity per unit stellar 
mass. The radii of the shells are larger than they would be in a more realistic model where 
the pressure of the ISM is included and where the effects of lower mechanical luminosity and 
leakage of bubble gas are included. These might be major effects and they both act in the 
same direction. If the shells at the best fit age are smaller, then we simply need more of them 
to build up the II2 mass sufficient to explain the observed CO flux. In particular, if the shells 
were to stall early at radii about half that in our model, then a model which includes this 
effect would require about eight times the cluster mass for the shell ensemble. This factor 
already accounts for much of the missing stellar luminosity /mass. Thus the stellar mass 
required is very sensitive to the adopted model. Additional simulations to test a stalled-shell 
scenario will be carried out in the near future. 

There are several other possible factors contributing to the shortfall in stellar luminosity 
represented by our model. Some SCs blow their shells out of the disk and are not detected, 
some or perhaps even most OB stars do not form in SSCs, and perhaps earlier generations 
of stars will augment the FIR luminosity to some degree. 

(4) Our model neglects the emission from the low density ambient ISM (rij^m = 30 cm~^), 
due to the lack of knowledge of the structure of this component in a starburst galaxy, and to 
the lack of direct observational data of this gas component that could be used to distinguish 
this gas and its physical state from the shell emitting gas in our models. If the ambient 
medium were uniform, as assumed in the model, it would produce no observable emission, 
since the density is too low to excite even the first excited rotational level. If, however, the 
ambient gas is assumed to be highly non-uniform, as is more likely the case, we can use the 
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total mass adopted for the sampled volume to estimate its CO emission by simply using the 
X-factor for the CO emitting gas in our own Galaxy. We find that, before the shells form, 
the ^^CO(l-O) emission from this ambient gas component would be about 54% of the total 
current emission within the central 1 kpc region. Note however that at high- J, there may be 
only very weak emission from this ISM component. Hence, the consequence of adding the 
emission of this lower density ambient gas component to the model would be to reduce the 
apparent excitation of the combined emission, especially the lower-J transitions. 

(5) A fixed microturbulent Doppler velocity {6vd = 1-5 km s^^) is used in our model, 
as a heating mechanism in the PDRs and as a broadening mechanism for the molecular line 
profile. However, the choice of turbulent velocity directly affects the computed CO line flux 
densities. In the optically thick case, the flux is directly proportional to the turbulent velocity, 
whereas in the optically thin case there is no dependence. Since the lower-J transitions are 
optically thick, more so than the higher level transitions, we anticipate that the use of a 
higher turbulent velocity would reduce the slope of the line SED (see Fig. [8]). 

(6) The assumption of the semi-infinite plane-parallel geometry in our PDR code is 
certainly a limitation, as the FUV intensity may be underestimated due to leakage of radi- 
ation from the region exterior to the cloud/shell, resulting in an increase in the local mean 
intensity at the edge of the slab. More advanced models of the shell geometry are simply 
beyond the scope of the PDR code at the time, and solving for the radiation field from both 
sides of the slab dramatically increases the computation time and would make the code too 
slow to run. 

(7) Other effects will invalidate our model for later stages of shell evolution. For example, 
after 30 Myr the largest radius of the shells in our model ensemble is about 678 pc. Thus, 
it will be merged with other shells, an effect which is not included in the model. It will also 
have extended beyond the scale height of the disk (300 pc along the minor axis), and be 
subjected to shear due to differential galactic rotation. However, these problems occur at 
ages older than our CO best fit model, and they should not significantly affect our best fit 
result. 

Overall, our analysis shows that the sub-mm/mm line emission reflects the recent star 
formation history in a starburst galaxy. The foregoing list of caveats and limitations ul- 
timately limit the precision with which one can obtain a realistic age for the starburst by 
the method described in this paper. Some of the effects described lead to an underestimate 
and some to an overestimate. Thus, to some extent, the effects are liable to cancellation. 
However, the one factor upon which the ages depend most strongly is stellar evolution, since 
the evolutionary state of the cluster governs the flux of FUV emission incident upon the 
shell, and this in turn has an important influence on the SED of the molecular line emission. 



-49- 



A consequence is that the age is unhkely to be profoundly affected by the effects hsted. 
This also means that there should be little surprise with the agreement with other meth- 
ods. However, the total cluster mass responsible is exceptionally sensitive to the model for 
the expansion of the shell, and consequently this quantity is less well determined than the 
age. In our model, it appears likely that this stellar mass in our cluster ensemble is severely 
underestimated. 



5. Conclusions 

This paper presents a first attempt at addressing the question of whether there is a 
signal in the FIR/sub-mm/mm molecular and atomic line data of the phase of a starburst. 
By comparing our evolving starburst models with available data of nearby starburst galax- 
ies, notably M 82, we show that it is possible to (1) successfully model the time-dependent 
FIR/sub-mm/mm line emission of molecular and atomic gas; (2) relate the observed molec- 
ular line properties of a starburst galaxy to its age, and hence to constrain the global star 
formation history; and (3) examine the possible relevance to the formation of the molecular 
rings in M 82. 

In essence, we have provided a complementary study to the previous work on estimating 
the age(s) of starburst in M 82 using quite different methods. In particular, the method is 
analogous to that of Efstathiou et al. (2000), which considered the observable effects of an 
evolving cluster on the IR emission from the surrounding expanding dust shell. Wc have 
also provided support for the hypothesis of molecular ring formation in the center of M 82. 

The main conclusions drawn from comparisons of our ESbM model with the observation 

are: 

1. There is good agreement between our supershell kinematic model and the observa- 
tions of the expanding supershell centered around the presumed supernova remnant SNR 
41.9 -I- 58 in M 82. The agreement supports the hypothesis that this supershell is created by 
strong winds from a young star cluster with a total mass of 2.2 x 10^ which formed at 
its center about 0.8 Myr ago, and the total mechanical energy needed for the creation of this 
supershell is about 1.5 x 10^^ ergs. This is the energy equivalent of the winds associated 
with ~ 1700 O stars (each with > 40 Mq). Our model also shows that there should be 
excess CO emission at high excitation transitions in this supershell. This is consistent with 
the provisional detection of such excess emission at ^^CO(6-5) in the region of this supershell 
seen after the surrounding disk emission is removed. Both agreements suggest that the set of 
evolving starburst models we have put forward in this paper can be used to interpret other 
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shells in M 82 or shells in other starburst galaxies. 

2. The age derived from our analysis of CO line SEDs and hne ratio diagrams using an 
instantaneous burst model is 5 - 6 Myr for the central 1 kpc region, with most of the CO 
emission arising from the central 680 pc region. The age derived from our atomic data is 
slightly older (10 Myr) for a larger area (i.e. 1,2 kpc) We suggest that these two ages may 
be a result of sampling different regions. It is unclear from our analysis whether they refer 
to the same period of star forming activity or to two spatially separated independent bursts. 
A more sophisticated model and more data are needed to clarify the picture. We do note 
however that our extended starburst model result also shows that a uniform star formation 
rate over the past 10 Myr can also produce a satisfactory model fit to the ^^CO data. Hence, 
the star formation in M 82 can be either viewed as instantaneous burst occurred 5-6 Myr 
ago, or this epoch could represent a characteristic time about which recent star formation 
history is centered. These burst ages derived from our models are similar to the results 
found in optical and infrared studies. These results lead us to conclude that the observed 
FIR/sub-mm/mm line spectra of a starburst galaxy can be successfully modeled in terms of 
the evolutionary scheme of an GMC/shell ensemble, and such studies can usefully constrain 
the age(s) or star formation history of a starburst galaxy. 

The starburst ages derived from our model are dependent on a great variety of assump- 
tions, e.g. the initial upper mass limit of the cluster spectrum. We find that the model 
can provide acceptable fits to the data only if the dominant initiating starburst clusters are 
massive, at least 5 x 10^ M©, corresponding to a GMC mass of 2 x 10^ M©. The uncer- 
tainty of the derived age is also affected by many other model assumptions, and the effects 
of varying these assumptions have not been examined. These include, for example, the effect 
of including the CO emission (especially low- J transitions) from the lower density ambient 
ISM, the effect of including a higher cosmic-ray ionization rate, and the effect of increasing 
the shell microturbulent velocity. Some of these effects would lead to an underestimate and 
some to an overestimate of the age, and hence to some extent, these effects would be ex- 
pected to cancel each other out. However, since the evolutionary state of the cluster governs 
the fiux of FUV radiation incident upon the shell, and this in turn has a pronounced effect 
on the SED of the molecular line emission, the stellar evolution is a crucial factor in con- 
straining the derived age. The starburst stellar mass and luminosity predicted by our models 
are significantly underestimated, based on a comparison with the observed FIR luminosity 
which is a factor of about twenty larger than our model value for the total luminosities of 
the clusters. Probable causes for this underestimate include (1) the neglect of the effects 
of the pressure exerted by the ambient gas, resulting in an overestimate of the shell radii; 
and (2) an overestimate of the supernova mechanical power which would also lead to an 
overestimate of the shell radii. Including these effects would allow more stellar luminosity 
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in the starburst for the mass of gas swept up in the shells required to match the CO data. 
Hence, the shortfall in our predicted stellar luminosity tends to support the wi dely held idea 



t hat bubbles / shells grow more slowly than the simple bubble theory predicts (IWeaver er al. 



19771 : iMcCrav fc Kafatoslll987l ). 



Our model also cannot provide a basis for incorporating higher density tracers (e.g. 
HCN and HCO"'"), Because the critical densities are, for example, approximately 10^ and 
10^ for the (1-0) and (4-3) transitions for both molecules. The HCO^ intensities may be 
affected by cosmic-ray ionization when it becomes the dominant heating source in the gas. 
Both HCN and HCO^ are associated with dense shells and their parent GMCs seen only at 
the earliest phase of the starburst evolution. But another principal source is required, most 
probably dense gas associated with the cores of potential star forming regions, not included 
in our model. If the dense cloud core component was included in our model, the effect of a 
higher cosmic-rate ionization rate (in M 82) should be taken into account. 

3. The results of the model analysis described above (item 2), also yield insights into 
the total gas content and its structure. For example, the total H2 gas mass ~ 2 - 3.4 x 10^ 
Mq, is consistent with that measured independently in the center of M 82. The inference is 
that the neutral ISM and possibly the molecular ring in the center of M 82 are largely the 
products of evolving shells. However, our interpretation concerning the ring formation is not 
unique, and the rings may also be created by Linblad resonance instabilities associated with 
the gravitational effects of the bar. 
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Tablc 3. Model parameters and variables 



Models 


Description 


Independent Variable: 


time or starburst age t 


Dependent Variables: 


shell radius Rs , expansion velocity Vs , number density Us , and thickness dg 




gas (kinetic) temperature Tgas, dust (thermal) temperature T^ust 




chemical abundances of different molecules and atoms in the shell 




number densities of coUisional partners H, e~, H+, P-H2, 0-H2, and He 


Fixed Parameters: 


CMC mass Mcmc- 3.16 x 10^ - 10^ Mq 




stellar mass m*: 0.1 - 120 Mq 




star formation efficiency (SFE) r] = 0.25 for Winds, r) = 1.0 for post-SN 




metallicity Z = 1.0 Zq 




gas-to-dust ratio = 100 




ambient ISM density of each shell nig^ (parent GMC at Winds, 30 cm~^ at post-SN) 




microturbulent velocity Svd = 1.5 km s~^ 


Fitting Pajrameters: 


total gas mass of the ensemble Mf^f^^i, burst age t 


Outputs: 


line profiles for each transition in each molecules and atoms 




integrated line intensity or flux 



